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We investigate the total asymmetric exclusion process by analyzing the dynamics 
of the shock. Within this approach we are able to calculate the fluctuations of the 
"^ . number of particles and density profiles not only in the stationary state but also 

j^ , in the transient regime. We find that the analytical predictions and the simulation 

H ' results are in excellent agreement. 
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1 Introduction 

The total asymmetric exclusion process (TASEP) is maybe the most simple lattice model 
of collective particle transport |jl|-p[. It can be interpreted, e.g., as a basic model of traffic 



^ [ flow m or as a model for biological transport g . In generic steady states, these processes 



are characterized by a non-vanishing mass transport which leads to an intrinsic non- 
equilibrium behavior. This is interesting from a theoretical point of view. For example, 
the stationary probabilities of the system are not given by an equilibrium ensemble. Of 
particular interest are systems with open boundary conditions because one can observe 
5^ I transitions between different bulk states induced by tuning the boundary conditions 0. 

Another reason for the importance of the TASEP model is the possibility to find an 
exact solution for the stationary state. This has been done by solving recursion relations 
for the stationary partition function [[^-[T0[| . Later a phenomenological approach, referred 
to as domain wall (DW) theory, was proposed |ll|] which was able to reproduce the exact 



results for the phase diagram, in the limit of large system sizes [|^. The interest of the 
latter theory is twofold: First, it leads to an intuitive understanding of the most important 
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features of the system and second it can be generalized to more complicated models which 



are out of scope for an exact treatment [|T^, ^ 



In our article we test the accuracy of the DW theory for the TASEP on finite open 
lattices. In particular we compare the stationary fluctuations of the number of particles, 
which are related to the integrated two-point correlation functions. Beyond that we 
also analyze the non-stationary behavior of TASEP, which cannot be treated with exact 
methods. Rather than investigating the relaxation spectrum of the process [|l5l, for which 



numerical evaluation is limited to very small system sizes (L < 10), we focus on the 
dynamics of directly observable quantities. We calculated e.g. the time dependent density 
profiles, which can be compared easily to simulations (which can be easily carried out on 
much larger system sizes). 

The outline of this paper is as follows. In section 2 we give a definition of the model 
and remind briefly of the most important features of the DW theory. Then we compare the 
predictions for the fluctuations of the number of particles with simulation results. Finally 
(section 4), we show how the system responds to a change of boundary conditions. 

2 The TASEP with open boundaries 

2.1 Definition of the model 

The TASEP is defined on a one-dimensional lattice of length L. Each lattice site {€) 
can be occupied by one particle (tj = 1) or be empty (tj = 0). In continuous time the 
dynamics of the particles is defined as follows: A pair of sites (i,i + 1) is chosen with 
probability dt, where dt denotes an infinitesimal time-step. If the site {i) is occupied and 
the site {i + 1) is empty one exchanges the positions of particle and hole. All other local 
configurations are unchanged. In case of open boundary conditions one additionally has 
to define the in- and output of particles. At a given time a particle can be introduced on 
the first site with probability adt if the first site is empty. If the last site of the chain is 
occupied the particle may escape from the chain with probability f3dt. 

The TASEP with continuous time dynamics can be implemented by a random sequen- 
tial update, i.e. one first chooses randomly a link between two sites and then performs the 
local update. A possible implementation of the open boundaries is to add two additional 
sites 0, L + 1 to the chain, which are occupied with probability a and 1 — (3 respectively, 
and act as particle reservoirs. By definition, one time step includes L + 1 link choices, 
and the corresponding local updates. 

For the TASEP it has been shown that the partition function, i.e. the sum of all 
stationary weights for a system of size L + 1, can be obtained recursively from the results 
for a system of size L [0,||. The recursion relations have been solved exactly P, e.g. by 
means of a matrix representation (MPA) |jTO|. Now we shall discuss how most of these 
results can be recovered using the DW theory. 



2.2 Domain wall theory 

The domain wall (DW) theory [0] gives a phenomenological description of the system 
dynamics. The basic idea of this picture is that, as long as the entrance (exit) capacity 
does not exceed the capacity of the chain, each particle reservoir enforces a domain (of 
constant density) in the bulk. At a given time both domains coexist in the chain. The 
coexistence of two domains in the chain implies the existence of a shock, i.e. a region 
where the two domains meet. In the domain wall theory one assumes that the shock is 
sharp, i.e. that it has a finite width W. For a large class of driven lattice gases this 
requirement is fulfilled and one has W '^ L already for moderate system sizes L. After 
the introduction of the DW theory for the TASEP it has also been successfully applied to 
models discrete in time and without particle hole symmetry |1^,|18|. Therefore the DW 
theory provides a quite general theoretical framework for models of particle transport. 

Below we now cite the ingredients of the DW theory which are of relevance for the 
remaining part of the article. In this theoretical frame, the dynamics of the shock de- 
termines the behavior of the system. A first characterization of the shock dynamics is 
possible by means of the lattice continuity equation 

^p(^,t)=j-_i(t)-j-(t), (1) 

where ji{t) denotes the local particle current at position i and time t and p{i, t) the density 
at the same site. By evaluating this in the continuum limit for the shock position one 
finds that the shock moves far from the boundaries with velocity 

V= ^+~^'- , (2) 

P+ - P- 

where p+{p-) and j-{j+) are respectively the densities and fluxes in the right (left) do- 
mains separated by the shock. 

The motion of the domain wall can be interpreted as a random walk with hopping 
rates 

D+ = , ^^ y D. = -^— (3) 

(P+ - P-) {P+ - P-) 

for a move to the right (left) ||TT|. In case of a blocked entrance or exit one can easily verify 
that the dynamics of the domain wall is described correctly. In order to generalize this 
observation, one uses that both particle reservoirs are independent. The interpretation of 
the shock dynamics as a random walk allows for an easy calculation of several quantities 
of interest. 

Now we specify the above quantities for the TASEP. An input rate a leads to a bulk 
density p- = a and an output rate /3 to p+ = 1 — /?. Inside each domain one obtains the 
same behavior as for the periodic system, e.g. one gets j = p(l — p) for the flow. 

In case of the TASEP one obtains still good results if one assumes that one can identify 
a single link as the position of the shock. The link is labeled i if localized between sites 



i and i + 1. Thus the wall location varies from to L for a lattice of L sites. On a finite 
system of size L the domain wall performs a random walk in a lattice with reflecting 
boundary conditions. Therefore the probability to find the domain wall at time t and 
position i can be evaluated from: 

^^^ = D+P{z-l,t) + D^P{z + l,t) 

-iD+ + D_)P{t,t), (4) 

for 1 < i < {L — 1). At the (reflecting) boundaries one has 

dP{0,t) 



dt 
dP{L,t) 

Jt 



Z}_P(l,t)-D+P(0,t) (5) 

D^P{L-l,t)-D_P{L,t), (6) 



where D^. = j^ — ^ and D_ = ^ — ^. 



The stationary solution of eq. (|H|) in the low density phase «</?,/?< 0.5 is given 
by P{i) = exp(— (L — x)/^)/A^, where the localization length is given by ^ = \og{D+/D_) 
and the normalization Af hy Af = (1 — exp(— (L + l)/^))/(l — exp(— 1/^)). In the high 
density phase one obtains analogous results by using the particle hole symmetry of the 
TASEP. By using choice (^) for the hopping rates, one recovers the exact result of the 
localization length ^. 

3 Particle number fluctuations in the domain wall 
picture 

The fluctuations of the number of particles A^ = ^^ Tj are defined as: 

{AN) {N^) - {Ny 



1 
L 



L 

(7) 



2j2{^.r,) + {N)-{N)' 



i<j 



where the brackets (...) denote an average over the stationary ensemble. Eq. (^ illus- 
trates the connection between the fluctuations and the two-point correlation functions. 
The fluctuations of the particle number, as well as the complete large deviation function. 



can be calculated using the exact stationary solution of the TASEP [|16| . We show now, 
as a sensitive test of the accuracy of the DW theory for finite systems, how the results of 
the domain wall theory compare to the simulation results. 



3.1 Particle number fluctuations and the domain wall theory 

We now proceed to calculate the relevant quantities in the framework of the DW theory. 
The first quantity one has to calculate is the average number of particles. In the low 
density phase the stationary probability to find the shock at link i is given by P{i) = 
(l/A/") exp(— (L— i)/(^). In the domain wall picture this means that one finds {L—i)p++ip_ 
particles in the systems. So by summing over all possible shock positions one gets 

(iV) = -lX]e-(^-)/^[(L-^)p+ + zp_] 

i=0 



+ 



p.L 

5e-i/e (1 - (L + 1) e-^/« + Le-(^+i)/«) 
(1 - e-i/?) (1 - e-(^+i)/5) ■ 



where S = p-^—p^ denotes the density difference between the high and low density domain 
which coexist. 

Next we have to calculate the two point correlations (tjTj) for i < j. This is done by 
averaging over the particle densities if the domain wall is located left from site i, between 
i and j, and right from site j, i.e. one finds 

fc=0 k=i 



k=j 



(9) 



P±\^(.~iL~k)/i_ 



Using eq. (0-^ we get the following result for the variance of the number of particles: 

{AN) 



P-(l-P-) 

6 (1 - 2p_)e-V^ fl-e^ _ _,/^ 
^ L l-e-(^+i)/« \l-e-y^ 

§2 g-2/?(i _ g-2L/€) _ 2e-(^+3)/?(l - e- 



L/^) 



(10) 



+ 



L (1- 6-1/5)2(1 _e-(L+i)/?)2 

1 + e-(^+i)/« + L 



52g-(L+l)/5. 



(1 - e-(^+i)/«)^ 



The corresponding results in the high density phase can be obtained by applying the 
particle hole symmetry of the TASEP. On the coexistence line {a = (3) the localization 
length diverges and the fluctuations are given by: 

^ = - (1 - 2p_) + (1 - p_) p_ + 6' ^—^. (11) 



3.2 Comparison with simulation results 
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Figure 1: Domain- wall prediction vs. simulation results of the particle number fluctua- 
tions {AN)/L. The symbols correspond to simulations results of the finite chain. The 
predictions of the domain wall picture are indicated by lines. We consider (3 = 0.25 for all 
simulations. The inset shows the relative deviations ((AA^) — (AN) dw) / (AN) between 
simulations and domain-wall predictions. For L > 300 the devations are smaller than the 
numerical precision of our simulations. 

Fig. m shows a comparison between the calculated values and simulation results for 
the particle number fluctuations. The results are for f3 = 0.25 and different values of 
a surrounding the first order transition at a = 0.25. Our simulations start in the low 
density phase, cross the coexistence line and end up in the high density phase. The 
agreement between DW-theory and simulation results is excellent for larger system sizes. 
For small system sizes one observes small deviations from the simulation results, which 
are presumably due to the finite width of the shock. Nevertheless, we stress the fact that 
the DW theory gives correct results (i.e. within the accuracy of the simulations) already 
for systems of the order of hundred sites. 



4 Non-stationary processes 

The results for the stationary particle number fluctuations presented in the previous 
section - for which an exact solution exists - have provided a test for the accuracy of the 
domain-wall picture. But in contrast to the exact stationary solution the domain wall 
picture also allows to characterize the dynamics of the system in the transient regime. 

4.1 Density profiles 

In order to check the correctness of the phenomenological picture we investigated the 
time-dependent density profile interplaying between two stationary states. We started 
from a stationary state in the low density regime on the line ao + /^o = 1- Indeed, on this 
line, the exact stationary solution has a simple product form, and the density profile is fiat 



T0| . This choice was made in order to keep the calculations simple, and to generate easily 
a large number of independent configurations. Suddenly, we change the output rate (3 
such as to (i) lie on the coexistence line (3 = a = a^ oi (ii) cross the transition line into the 
high density phase so that (3 < a = a^. We take as the origin of time the moment when 
this change occurs. A large number of independent simulations are performed in parallel. 
At regular time intervals, an ensemble average of the density profile over the independent 
simulations is performed. After some time, a new stationary state is reached. 

In our simulations, we chose as initial condition oq = 0.25, /3o = 0.75. At t = 0, 
a = ao stays fixed, while /3 is varied to (i) (3 = 0.25, or (ii) /3 = 0.2, 0.1, or 0.02. This 
change is instantaneous. The density profile is averaged every 100 time steps over 100,000 
independent simulations, i.e. we took 100,000 arbitrary configurations from the initial 
stationary ensemble and performed a simulation run for each initial configuration. 

In the remaining part of this section we want to explain how this setup can be described 
by means of the domain wall theory. 

The chosen set of parameters ensures the absence of a shock in the initial condition. 
By changing the right reservoir we introduced a shock at position L and time t = 0. In the 
course of time the motion of the wall then follows a random walk described by equations 



In case (i) the random walk is symmetric, because a = P and thus D^ = D_. An 
exact analytic solution for a random walk between two reflecting walls, starting at r=0. 



is known for all times 17 



P(^,t) = _l-|l + |^exp(-2Dt|l-cos[-^]|) 



/ inn \ /'(^ + 1)7™ 

cos + cos 



L + lJ \ L + 1 



(12) 



At long times, the wall has a uniform probability to be located anywhere in the system. 
In case (ii), D^ and D_ are different. The exact analytical solution of the asymmetric 
random walk with reflecting boundaries is also known |[l^,|T^, but for convenience we 



treated directly the discretized diffusion equation 



P{i,t + dt) = D+dtP{i-l,t) + D_dtP{i + l,t) 
+ [l-{D+ + D_)dt]P{i,t), 



(13) 



with equivalent expressions for the boundaries. We took dt small enough so that the 
results do not depend on the chosen value. 

Once the distribution P{i,t) is known, the density profile can be computed for any 
time 



p{z,t)= K^P(n,t) p_+ J2 ^(^'i) P+ 



(14) 



,n=0 



\.n=i+l 



where p_ and p+ are the effective densities of the left and right reservoirs, i.e. p- = a 
and p+ = 1 - /? for the TASEP. 

a and (3: (0.25,0.75)->(0.25,0.25) 




Figure 2: Domain- wall prediction (dashed lines) vs. simulation results (solid lines) for 
the time dependent density profile when /? is suddenly changed at t = from /3 = 0.75 
to 0.25. The input rate is fixed a = 0.25 and the system size is L = 100. Each profile 
was averaged over 100000 independent simulations. The dot-dashed line indicates the 
asymptotic linear density profile. 

The comparison with the simulation results (see figs. ^ and |^) shows an excellent 
agreement. In case (i), the domain wall undergoes a diffusive motion and the density 
profile modification scales as \/t. To make the picture more clear, beyond t = 500, we 
show only the profiles every 500 time steps. In the long time limit, the domain wall is 
delocalized over the whole system and the density profile is linear. 
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a and (3: (0.25,0.75)->(0.25,0.1) 
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Figure 3: Domain- wall prediction (dashed lines) vs. simulation results (solid lines) for 
the time dependent density profile when /3 is suddenly changed at t = from (3 = 0.75 
to 0.1. The input rate is fixed a = 0.25. Profiles are plotted every 100 time steps. Each 
profile was averaged over 100000 independent simulations. The system size is L = 100. 

In case (ii), the wall has a non vanishing drift velocity, and the shift of the density 
profile is linear in time, at least as long as boundaries are far enough. This is illustrated 
in fig. ^ where the first density profiles (for 100 < t < 500) are translated by a multiple 
of 100 V. The averaged wall velocity V is given by eq. (^. Of course, as the wall position 
distribution evolves from a Dirac function into a larger and larger gaussian, the density 
profile flattens as time increases. At long times, the wall is located near the left boundary, 
and the density profile is fiat, apart from an exponential boundary layer at the entrance 
of the system (fig. |) . 

For small systems, the agreement between domain wall theory and simulations is not 
so good, as expected. However, even for a system of size L = 10, the domain wall results 
still indicate quite well the non- stationary behavior, as long as the wall is far enough from 
the left boundary (fig ^). Surprinsingly, the agreement is much better at the beginning, 
when the wall is still near the right boundary. It could be linked to the fact that in the 
initial condition, as the density profile is fiat, the domain wall approximation (with a wall 
on the right of the system) is exact. 



4.2 Fluctuations in the particle number 

So far we focused on density profiles. Now we shall show that the domain wall picture 
allows also to predict the fluctuations in the total number of particles in a dynamical 
regime. From the calculation in the previous section, the time dependent probability 



a and P: (0.25,0.75)->(0.25,0.1) 




Figure 4: Same figure as H, but with the density profiles translated by a multiple of 100 
V. 



distribution P{i,t) is known. The fiuctuations are evaluated from formulas similar to 
), but now we replace the stationary distribution e~^^~^''^ by P{i,t). This yields 



(iV)(t) = Lp^ + (p+ - p_) J2iL - k)P{k, t) 



(15) 



fc=0 



and 



2 *— 1 i~i 

fc=0 k=i 



.2 i 
k=j 



+ '^Ep(Kt). 



(16) 



Using these expressions with the definition (|^), the variance of the number of particles can 
be computed numerically for any time, and compared with the direct simulation results. 
Fig. P presents such a comparison in case (ii), i.e. for a final value /3 = 0.1 Time has been 
divided by the average time needed for the wall to cross the system, i.e. L/V where V is 
given by equation (^. The agreement is excellent for all sizes, though for L = 50, some 
small finite size effects are visible. 

In the early linear stage, the wall has not yet reached the left boundary, and the curves 
superimpose for all sizes. When the wall arrives near the left boundary, the width of the 
probability distribution P{i,t) scales as vL. Thus the typical time during which domain 
walls reach the left boundary in the different realizations of the system also scales as vL. 
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Figure 5: Domain- wall prediction (dashed lines) vs. simulation results (solid lines) for 
the time dependent density profile when /3 is suddenly changed at t = from /3 = 0.75 to 
0.1. The input rate is fixed a = 0.25. Profiles are plotted at time 0, 5, 10, 15, 20, 25, 28, 
31, 34, 37, 40, and 500. Each profile was averaged over 100000 independent simulations. 
The system size is L = 10. 

As time has been rescaled by L in fig. |^, this transient separating the linear growth from 
the final stationary state tends to a sharp jump as L becomes large. 

In the final stationary state, the fluctuations within the high density phase dominate. 
They can be estimated from the mean field expression p+(l — p+) = (1 — /?)/? = 0.09. 



5 Conclusion 

Our results have shown that many features of the TASEP with open boundary conditions 
can be understood by analyzing the domain wall dynamics. For the stationary state, 
where the DW theory is known to give exact results in the limit of large L — ^ oo |T^ , we 
calculated the fluctuations of the particle number for chains of various sizes. As being 
the integrated two-point correlation functions the fluctuations of the particle number are 
a sensitive test for the accuracy of the DW theory. We found that the predictions of the 
DW theory are numerically indistinguishable from the simulation results for L > 100. 

But beyond that we also found an excellent agreement in the transient regime. This is 
of particular interest, because the transient regime is out of scope of exact solution. The 
perfect agreement (up to the precision of our measurements) between simulation results 
and DW theory shows that the domain wall motion determines exclusively the relaxation 
processes, i.e. the coupling to the particle reservoirs is instantaneous. For future work it 
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Figure 6: Fluctuations in the total number of particles as a function of time (rescaled 
by L/V). The simulations (symbols) are compared with the domain wall estimate (black 
continuous lines) with excellent agreement. The dashed line indicates the value p+(l — 
p+) = 0.09. 

will be interesting to see whether other relaxation mechanisms are of importance if the 
domains are not of a simple product measure. 

In several other studies it has been shown that the DW approach correctly predicts the 
phase diagram for a large number of models and update schemes [0, 0, 0. In general 
cases, where an exact solution of the process is not possible, one has to establish the 
coupling to the chain and the flow density relations numerically. However, also in this 
case the domain wall interpretation is useful: First of all the parameter space is reduced 
drastically (e.g. the bulk density corresponding to a particular input prescription can be 
obtained from a single simulation run with free right boundary), and second it allows for 
a better characterization of the physical behavior [l^. Let us illustrate the latter point 
by the example of the Nagel-Schreckenberg model for traffic flow with maximal velocities 
Vmax > 1 mi- In this case one has several possibilities to implement the in- and output 
of the particles. This choice can produce even qualitatively different results for the phase 
diagram. E.g. it is necessary to apply a reservoir that can achieve the full capacity of 
the chain for a given model, in order to observe the maximum current phase (see |2^ 
for a counter-example). While for the TASEP the natural input of particles on the first 
site of the chain fulfills this prerequisite, this is more subtle for particles which can move 
with higher velocities [jl4|- The difficulties in finding a suitable particle reservoir can be 



circumvented using the basic idea of the domain wall theory for the description of the 
particle reservoir. Instead of describing the capacity of the reservoirs in terms of in- and 
output probabilities we recommend to describe a particle reservoir by its effective density. 
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i.e. by the bulk density which corresponds to a specific parameter combination, and to 
calculate the phase diagram only thereafter. 

Summarizing, we have shown that the DW theory is able to reproduce the behavior 
of the TASEP with open boundary conditions to a large extend, and we have used it to 
obtain new results for non- stationary fiows. We expect that our results can be generalized 
to a wide class of models for particle transport. 
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